Spin 1/2 Fermions in the Unitary Regime: A Superfluid of a New Type 
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We have studied, in a fully non-perturbative calculation, a dilute system of spin 1/2 interacting 
fermions, characterized by an infinite scattering length at finite temperatures. Various thermo- 
dynamic properties and the condensate fraction were calculated and we have also determined the 
critical temperature for the superfluid-normal phase transition in this regime. The thermodynamic 
behavior appears as a rather surprising and unexpected melange of fermionic and bosonic features. 
The thermal response of a spin 1/2 fermion at the BCS-BEC crossover should be classified as that 
of a new type of superfluid. 
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PACS numbers: 03.75.Ss 

The unitary regime is commonly referred to as the situ- 
ation in which the scattering length a greatly exceeds the 
average inter-particle separation, thus n\a\ 3 ^> 1, where 
n is the particle number density 0,0- It is widely ac- 
cepted by theorists that at T — these systems are su- 
perfluid and that in the unitary regime the coherence 
length is comparable in magnitude with the average in- 
terparticle separation. At T = this problem has been 
considered by a number of authors 3] and the most ac- 
curate results so far have been reported in Refs. El- 
la 2002 it was shown experimentally that such systems 
are (meta)stable, and they have been studied extensively 
experimentally ever since 

The typical theoretical treatment of such systems is 
based on the idea put forward by Eagles, Leggett and 
others |{J, and used subsequently by most authors [ToL 
111) . The form of the many-body wave function is as in 
the weak coupling BCS limit and is used for all values of 
the scattering length a. The particle number projected 
BCS wave function has the functional form 

*(ri,r 2 ,r 3 ,r 4 , ...) oc A[^>(r 12 )(j)(r 34 )...], 

where odd subscripts refer to spin-up particles and 
even subscripts to spin-down particles, A is the anti- 
symmetrization operator, r±2 = |ri — r 2 | and 4>(r) is ei- 
ther the Cooper pair wave function in the BCS limit, 
or the two-bound state wave function in the BEC limit. 
The main difficulty with this approach becomes evident 
when one tries to use this kind of wave function in the 
unitary regime, where n\a\ 3 3> 1. In the extreme BEC 
limit, this wave function describes a state with all bosons 
(dimers) at rest, in the condensed state. The fraction of 
non-condensed bosons (dimers) is known to be small then 
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where rid = n/2 and add — 0.6a is the dimer-dimer scat- 
tering length 0, • When one approaches the unitary 



regime, the fraction of non-condensed bosons becomes of 
order one fl3j . which resembles qualitatively the situa- 
tion in superfluid 4 He, and then a meanfield description 
(with or without fluctuations) becomes questionable. 

In order to calculate the thermal properties of a sys- 
tem of fermions in the unitary regime, we have placed 
them on a 3D-spatial lattice and used a path integral 
representation of the partition function. We start from 



Z(fl, n) = Tr [] cx P [-t(H - fxN)} 
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0(j3,fi) = 



^Tr|oncx P [-^- M A)]|, (2) 



where j3 = l/T = N t t and O is a quantity of interest. T 
stands for the temperature and [i for the chemical poten- 
tial, and H and N are the Hamiltonian and the particle 
number operators respectively. 

Since the system under consideration is dilute, we shall 
use a zero-range two-body interaction with a cut-off in 
momentum. Specifically, V(ri — r 2 ) = — g8(ri — r 2 ), 
with the additional prescription that all two-body ma- 
trix elements of this interaction vanish, if the relative 
momentum of the two particles exceeds a given cut-off 
momentum hk c . The renormalized coupling strength of 
this interaction is given by the following prescription 
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As we ultimately place the fermions on a spatial 3D- 
lattice, an implicit cut-off momentum is introduced by 
the lattice spacing I. If we were to allow the particles to 
move without restriction in all three spatial directions on 
such a 3D-lattice, we would be able to solve exactly the 
quantum mechanical problem with the restriction only on 
the particle momenta |14j . We impose periodic boundary 
conditions and consider the many fermion system in a 
cubic box of side L = NJ. It would be desirable to 
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have L exceed significantly the coherence length. At T = 
this condition is easily satisfied for a many fermion 
system interacting with a large scattering length, when 
n\a\ 3 ^> 1. At temperatures \T— T c \ -C T c this is not true 
anymore, as the coherence length diverges when T — > 
T c . The phase transition on a finite lattice is rounded 
and it will not show the expected singular continuum 
behavior. Since the momentum space on a 3D-latticc 
has the shape of a cube, while typically in field theoretical 
models with momentum cut-off the shape is spherical, we 
have included in calculations only 3D-momenta satisfying 
the condition k < k c < ir/l, in order to facilitate the 
analysis. 

The next step is to generate a sufficiently accurate rep- 
resentation of the propagator exp[— t(H — /J,N)] as 
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exp[-r(i? - fiN)] 
t(K- [j,N) 



exp(— tV) exp 



t(K - fiN) 



where K is the kinetic energy operator. The inter- 
action becomes a simple Hubbard attractive potential 
V = —gJ2i "4 where i labels the 3D-lattice sites 

and n^iii) are the number densities for the two spin 
states at a given spatial site i. The action of each fac- 
tor is evaluated either in coordinate or momentum space 
respectively, and the Fast Fourier Transform is used to 
connect these representations. The kinetic energy has the 
correct dispersion for momenta smaller than the cut-off 
momentum Hk c , namely £k = h 2 k 2 /2m. We have used 
a discrete Hubbard-Stratonovich representation of this 
interaction energy, similar to Ref. |l5j | 

exp[ffrri T (i)n|(i)] = 

\ J2 [l + MU)fh(W+MU)ni(i)]- 

<r(i) = ±l 



Here A = \J exp(<7r) — 1 and j is the label for the cor- 
responding imaginary time step. The partition function 
can then be expressed in the form 

Z(fi,n)=Tc J l[Va(i,j)U({a}), 

U({a})=T T eM-r[HW})-tA} 

where T T stands for a time-ordering operator. h({a}) 
is a one-body Hamiltonian, which has exactly the same 
form for both spin-up and spin-down states, and the mea- 
sure Y\ij "D <r(i, j) is assumed to contain the appropriate 
normalization factors. The expectation value of any op- 
erator takes the form 
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TrW({a}) 



(3) 



The trace can now be evaluated over all possible Slater 
determinants with various particle numbers and various 



expectation values acquire very simple forms 0, Il7| . 
There is no fermion sign problem in this case, since 



Tr U{{a}) = {det[l + U{{a})}} 2 > 0, 



(4) 



and where the determinant is computed in the spin-up (or 
spin-down, which is identical to spin-up) single-particle 
Hilbert space. The many-fermion problem has thus been 
reduced to a typical auxiliary field Quantum Monte Carlo 
problem, to which the standard Metropolis algorithm can 
be applied, using Eq. (0} as a measure. At each Monte 
Carlo step we have randomly changed the signs of a frac- 
tion of the cr-fields at all spatio-temporal lattice sites. We 
have increased or decreased the fraction of sites where 
the cr-fields were updated, so as to maintain a running 
average acceptance ratio (over the latest hundred Monte 
Carlo steps) between 0.4 and 0.6. In our simulations we 
evolved in imaginary time single-particle wave functions 
(plane waves) with momenta hk < hk c and calculated 
the measure Eq. @. In order to avoid numerical insta- 
bilities at low temperatures the Singular Value Decom- 
position technique was used |l6| . The expectation values 
in Eq.(PJ) were computed using the one-body density ma- 
trices 



™T,l( x ;y)= ^(x) 
ki,k 2 </c c 



U({*}) 



l+U{{a}) 



ki,k 2 



where ^k(x) = exp(ik • x)/L 3 / 2 . 

M. Wingate 01 j using the formalism of Ref. [l^. has 
estimated the critical temperature to be T c « 0.05 £f, 
but for a value of the scattering length a that was not 
determined precisely. The similar treatment in Ref. ji^l 
has in our opinion large discretization errors. In both 
approaches the choice of the kinetic energy as a simple 
hopping term can also lead to significant systematic ef- 
fects. 

The results of our simulations for lattices ranging from 
6 3 x 1361 and 8 3 x 1732 (at low T) to 6 3 x 300 and 8 3 x 257 
(at high T) and for 2 • • • 20 x 10 5 Monte Carlo samples (af- 
ter thermalization) are shown in Fig. ^ The imaginary 
time step was chosen as r = min(m? 2 /157r 2 S 2 , In 2/10g). 
We have estimated that the Monte Carlo correlation 
length is approximately 150 Metropolis steps at T w 
0.2ef- Consequently, the statistical errors are of the or- 
der of the size of the symbols in Fig. ^ The chemical 
potential was chosen in such a way as to have a total of 
about 20 particles for the 6 3 lattice and about 55 parti- 
cles for the 8 3 lattice. In all runs the single-particle oc- 
cupation probabilities for the highest energy states were 
significantly below a percent at all temperatures. 

At T c = 0.23(2) e F the behavior of E(T) changes. 
Here Ef is the Fermi energy of the free Fermi gas with 
the same number density n — N/L 3 . Mcanficld plus 
fluctuations estimates put T c at values slightly above 
the condensation temperature in the BEC limit, namely 
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FIG. 1: The total energy E(T) is shown with open circles for 
a 8 3 -lattice and with triangles for a 6 3 -lattice, and the chemi- 
cal potential fi(T) with squares for the case of a 8 3 lattice. The 
combined Bogoliubov- Anderson phonon and fermion quasi- 
particle contributions E p h+ qp (T) (Eq. |SJ is shown as a 
dashed line. The solid line is E Fg {T) - 0£e F N(l - f„), 
where EF g (T) is the energy of a free Fermi gas. In the up- 
per left inset we show the entropy per particle S(T)/N = 
[5E(T) /3 - n{T)N]/NT with circles for 8 3 and squares for 6 3 
lattices respectively, and with a solid line the entropy of a free 
Fermi gas with a slight vertical offset. In the lower right inset 
we plot the condensate fraction a(T) as defined in Ref. fl3|l . 
with circles the 10 3 -lattice results, with squares the 8 3 -lattice 
results and with triangles the 6 3 -lattice results, and the solid 
curve is a(T) = a(0)[l - (T/T c ) 3/2 ]. 



Tbec — 0.218 Ef 1 1 j - If a Fermi gas exactly at reso- 
nance behaves as a BCS superfluid, then its critical tem- 
perature would be T c 0.277 sp, when including the 
correction due to Ref. ; 22j. A meanfield plus fluctua- 
tions approach predicts T c sw 0.26 Sf |2lj . 

At very low temperatures one can expect that only 
two types of elementary excitations exist, the boson- 
like Bogoliubov-Anderson phonons and the fermion-likc 
gapped Bogoliubov quasi-particles. One can estimate 
their contribution to the total energy E(T) by assum- 
ing that at T = the system is a Fermi superfluid with 
a ground state energy and a pairing gap determined in 
Ref. H: 



Eph+qp{T) 




A 



2\ 7/3 
- £Fexp 
e / \2kpa 



(5) 



(6) 



determined in Ref. |4j to be very close to the weak cou- 
pling prediction of Gorkov and Melik-Barkhudarov [22^ . 
and £ s » 0.44 and Bp — h 2 kp/2m and n — k F /3ir 2 re- 
spectively. The sum of these contributions is plotted in 
Fig. as a dashed line. Numerically, both these contri- 
butions are comparable in magnitude over most of the 
temperature interval (0,T C ). One should not take seri- 
ously the apparent agreement though, as these expres- 
sions are only approximate formulas for T <C T c . It is 
notable that at temperatures in the vicinity of T c both 
fermionic and bosonic elementary excitations seem to be 
equally important, unlike the BCS or BEC limits, where 
only one type of excitations is relevant. At T > T c the 
system is expected to become normal. The fact that its 
specific heat is essentially that of a normal Fermi liquid 
Ep g {T) is however somewhat of a surprise, as one would 
expect the presence of a large fraction of non-condensed 
unbroken pairs. The caloric curve for a free Fermi gas 
Ep g {T) was offset vertically, so as to agree approximately 
at T = with the estimate of the energy of a Fermi gas 
in the unitary regime in the normal phase. The value 
can be estimated using the (approximate) condensation 
energy w -3A 2 N/8e F - 

One cannot fail to notice the behavior of the chemical 
potential (J-(T), which is essentially constant for T < T c , 
and decreasing with T, as expected, for T > T c . Apart 
from the natural vertical offset, this behavior is very sim- 
ilar to the behavior of the chemical potential for an ideal 
Bose gas undergoing condensation and it has some un- 
expected consequences. The T-dependence of the energy 
of a spin 1/2 fermion system at unitarity can be repre- 
sented by introducing the universal function £(x) (with 

m = e.) as 



E(T) = N^ept (— 
5 \e F 



which together with \i = const at T < T c implies 



(7) 



e (£)= 6 + c(£) , where n = \. (8) 

This temperature dependence is characteristic of an ideal 
(sic!) Bose condensed gas, even though the system is also 
superfluid at the same time. From our simulations we 
infer that the value of the exponent cannot differ from 
n = 5/2 by more than about 10%, and that values either 
n < 2 (n — 2 would be expected for a normal Fermi 
system) or n > 3 (n — 4 would be expected at T <C T c 
for a fermion superfluid) are inconsistent with our data. 
Our results are consistent with an effective boson mass 
to* 3to in this temperature interval (determined from 

COC TO* 3 / 2 ). 

One can also show that the entropy is given by 



where A is the approximate value of pairing gap at T = 



o rT/sF 

S(T) = -N dy 



(9) 
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where the prime indicates a derivative with respect to 
the argument. Thus S(T) oc T 3 / 2 for T < T c , which is 
intermediate between the behavior of a Fermi (oc T) and 
a Bose/phonon (oc T 3 ) systems. Since S(T) in either the 
BEC or BCS limits can be easily determined |2^|, one 
can use the entropy S(T) calculated here, see Eq. © 
and the upper left inset in Fig. (JTJ, to construct an ab- 
solute temperature scale in the unitary regime, using an 
adiabatic tunning of the scattering length, and extending 
thus the ideas of Ref. to the unitary regime. 

Apart from various themodynamic potentials, we have 
computed the temperature dependence of the condensate 
fraction ot(T), as defined in Ref. 0] and evaluated for a 
r = L/2 pair separation. The condensate fraction defines 
the off diag onal long range order of the two-body density 
matrix (24) . In complete agreement with the behavior of 
the thermodynamic potentials, the temperature depen- 
dence of the condensate fraction a(T) is consistent with 
T c = 0.23(2)ei?. The functional form of a(T) is, surpris- 
ingly, similar to that of an ideal Bose gas, see right inset 
in Fig. (JTJ. 

The value of T c determined here cannot be compared 
with the recent experimental result from Duke University 
H . The presence of a trap can change significantly the 
thermodynamic properties, and, in particular, the surface 
modes can play an unexpectedly large role, see Ref. |2fij . 

In conclusion, we have performed a fully non- 
perturbative calculation of the energy of a spin 1/2 
fermion system in the unitary regime, by placing the par- 
ticles on a judiciously chosen spatial 3D-lattice, in a path 
integral formulation of the many-body problem. We have 
determined the critical temperature of the superfluid- 
normal phase transition as T c = 0.23(2) ef- The ther- 
modynamic behavior appears as a rather surprising and 
unexpected melange of fermionic and bosonic features. 
The thermal response of a fermion system at the BCS- 
BEC crossover suggests that such a system should be 
considered a new type of superfluid. 

Discussions, and help with various computing issues, 
with Y. Alhassid, G.F. Bertsch, D.B. Kaplan, M. Prange, 
J.J. Rehr, B. Sabbey, D.T. Son, B. Spivak, M.B. Wingate 
and Shiwci Zhang are gratefully acknowledged. Support 
from the Department of Energy under grant DE-FG03- 
97ER41014, from the Polish Committee for Scientific Re- 
search (KBN) under Contract No. 1 P03B 059 27 and 
the use of computers at the Interdisciplinary Centre for 
Mathematical and Computational Modelling (ICM) at 
Warsaw University are appreciated. 

Note added. Burovski et al |2fij performed a finite size 
analysis and claim that T c /sf = 0.152(7). While we 
believe that the system sizes used are too small to evi- 
dence the critical power law behavior near T c , we notice 
that their results are in remarkable agreement with ours 
and, combined with T = data of Refs. [4,6], show that 
C(T) = AE(T)/AT is largest for T/e F G (0.15,0.3), 
which is roughly consistent with our findings. 
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